Supersolids versus phase separation in two-dimensional lattice bosons 
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We study the nature of the ground state of the strongly-coupled two dimensional extended boson 
Hubbard model on a square lattice. We demonstrate that strong but finite on-site interaction U along 
with a comparable nearest-neighbor repulsion V result in a thermodynamically stable supersolid 
ground state just above half-filling, and that the checker-board crystal is unstable for smaller V, 
and for any V just below half-filling. The interplay between these two interaction energies results 
in a rich phase diagram which is studied in detail using quantum Monte Carlo methods. 



The detection of possible supersolid (SS) state in re- 
cent experiments on solid ''He by Kim and Chan |l) has 
led to renewed interest |^ in a problem that has long 
1^ intrigued physicists: Can a supersolid phase — with si- 
multaneous diagonal (solid) and off-diagonal (superfluid) 
long-range order — exist in a bosonic system? While the 
the issue remains controversial ^ ^ in a translation- 
ally invariant system despite almost fifty years of theo- 
retical research, the situation in lattice models is clearer. 

Theoretical studies § H |] of various lattice boson 
models (which can nowadays be implemented using cold 
bosonic atoms on optical lattices ||l^), appeared to con- 
firm that here supersolid ground states can indeed exist, 
particularly when doped away from half-filling. Studies 
of the closely related quantum phase model found su- 
persolid order in the ground state even at half-filling 
However, as was pointed out recently the stability of the 
supersolid against phase separation had not been inves- 
tigated 0. Indeed, for hard-core bosons on a square 
lattice, the most widely discussed supersolid pattern — 
with (tt, tt) diagonal order — is thermodynamically unsta- 
ble and phase separates into a pure (tt, tt) solid and a 
superfluid (SF) for all values of interaction strengths. A 
striped supersolid phase — with (0, tt) ordering — is stabi- 
lized by a finite next-nearest-neighbor (nnn) interaction. 

In this work we analyze stability of crystalline and 
supersolid orders of lattice bosons. We present exact 
strong-coupling arguments showing under which condi- 
tions checkerboard supersolids are unstable, and how 
they can be stabilized with large but finite on-site and 
nearest- neighbor (nn) energies U and V. We support 
these arguments by quantum Monte Carlo simulations of 
a two-dimensional (2D) extended Bose Hubbard model 
demonstrating that the supersolid (SS) phase is stabilized 
for densities p > 1/2 and sufficiently large V (Fig. ^. 

Specifically, we study the extended Bose-Hubbard 
model (EBHM) on a d-dimensional hypercubic lattice 
with on-site (U) and nn (V) interactions. 
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FIG. 1: The ground state phase diagram of the 2D extended 
Bose-Hubbard model (|l]) in the V ~ p plane for U/t = 20 and 
densities p < 1, showing superfluid (SF) phases, checkerboard 
solids formed by single bosons (CDW I) and pairs of bosons 
(CDW II), a Mott-insulating phase (MI), phase separation 
(PS) and finally a supersolid phase (SS). 

where al{ai) creates (annihilates) a boson at site i with 
the occupation number rii = a^a;, t is hopping, fj, is the 
chemical potential, and (i, j) runs over all nn pairs. 

In the zero-hopping limit, t = 0, the non-negative po- 
tential energy {U, V > 0) is minimized at half-filling, 
p = 1/2, by the crystal state with only one sublattice 
occupied [checkerboard pattern with (tt, tt) modulation 
in 2D]. This state is gapped; it remains stable in the 
presence of a small hopping, t ^ U,V, with a kinetic en- 
ergy gain AE « —zt^ /[{z~\)V] per boson, where z = 2d 
is the coordination number. 

Introducing holes only costs chemical potential fi but 
no potential energy; the kinetic energy gain is somewhat 
increased but remains quadratic in t for isolated holes. 
However, the kinetic energy gain becomes linear in f if a 
number of holes encircle a region of a crystal [Fig. ^(a)]. 
The energy gain is maximized at AE ct,l<c<2 
per hole for a planar [linear in 2D, see Fig. ||(b)] domain 
wall doped with one hole per two sites. As a result, for a 
large system with N = L'^ sites, the crystalline order is 
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FIG. 2: The p — 1/2 checkerboard crystal doped with holes, 
(a) Four holes encircle a boson which can hop between the five 
degenerate sites, (b) Domain wall doped with holes; bosons 
can hop freely across the dashed line. 

destroyed by introduction of a small density p ^ of 
holes. This instability of the p = 1/2 crystal to domain 
wall formation upon hole doping excludes the possibility 
of a SS phase. In practice, on the isotropic square lat- 
tice, the instability develops further, leading to a phase 
separation between the commensurate crystal at p = 1/2 
and a uniform superfluid with p < 1/2 (Fig. |l|). 

Doping of the p = 1/2 crystal with additional bosons 
works differently depending on the relation between V 
and U. The energy cost to place a boson at an empty (oc- 
cupied) site is Eo = zV — p {Ei = U — p). Respectively, 
for U > zV, the additional bosons fill empty sites and 
mask the checkerboard modulation; for U — zV ^ i > 
the situation is precisely particle-hole conjugate to hole 
doping. The kinetic energy is again minimized at planar 
domain walls which destabilize the checkerboard crystal 
order. In particular, in the hard-core limit \] — > oo, the 
crystalline order is always unstable for p^ 1/2. 

With zV > U, however, the bosons (unlike holes) can 
be placed on either an occupied or unoccupied site. The 
total energy of a single boson delocalized between the 
two sublattices is S = E^i + A - {2zH'^ + A^)!/^, where 
A = (zV - U)/2. Clearly, for sufficiently small A - 
the kinetic energy E ~ Ei is again linear in t and large, 
which prevents the domain wall formation. As a result, 
these doped particles will form a superfluid on top of the 
charge ordered background and hence a supersolid. Two 
bosons experience both on-site and nn repulsion (2U and 
V respectively). Therefore, at sufficiently small densities 
the condensate should remain stable, which completes 
the formal argument for the supersolid existence. 

Similarly, at unit filling, p = 1, the ground state is a 
Mott insulator with one boson per site for U > zV, and 
an ordered solid with two bosons on every other site for 
U < zV. In the former case, additional holes (particles) 
move along the uniform background with the hopping 
integral t (2t) and experience both nn and on-site re- 
pulsion (infinite in the case of holes). They condense 
on top of the uniform background forming a superfluid. 
However, for zV > U, the doped particles move on a 
checkerboard background with the effective hopping, e.g., 

= 2t'^/[zV -U+{z- 2)V] for holes. The resulting ki- 
netic energy gain is only quadratic in t and can be super- 
seded if the holes come together into a supersolid phase 
with p > 1/2. Overall, this leads to a thermodynamical 
instability of the hole-doped checkerboard solid formed 
by pairs of bosons at p = 1: the system can minimize its 



energy by phase separation. Note that here phase separa- 
tion is not between a superfluid and a solid but between 
a supersolid and a solid. The solid order is not desta- 
bilized at this first order phase transition, but just the 
"Bose-Einstein condensation" transition of holes doped 
into the solid becomes first order. 

We next perform quantum Monte Carlo simulations 
to corroborate these arguments and to show the phase 
diagram and the existence of a supersolid phase for the 
EBHM in the low-density region p < 1 . 

We have used loop-operator updates in a stochastic 
series expansion (SSE) quantum Monte Carlo (QMC) 
method ^ to study the EBHM (|l|) in the strong cou- 
pling regime {U,V t) and for < p < 1. In the present 
study, simulations have been carried out in square ge- 
ometry, N = L X L, with L = 6, . . . , 16. Ground state 
properties have been obtained by taking sufficiently large 
values of the inverse temperature (3, where (3 = 2L turned 
out to be sufficient. 

To characterize different phases, we have studied the 
static staggered [Q = (tt, tt)] structure factor. 
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which measures the diagonal long range order (checker- 
board solid) in the system, and the superfluid density 
Ps, measured from the winding numbers of the bosonic 
world lines [Wx and Wy) in the x- and y- directions as 
Ps ~ {W]: + Wy)/2j3m, where m — 2/t \s the effective 
mass of the bosons. A checkerboard solid ground state 
at p = 0.5 is marked by a diverging S{tt, tt) and vanish- 
ing Ps, whereas a pure superffuid phase has S'(7r,7r) = 
and Ps > 0. A supersolid phase, on the other hand, is 
characterized by a diverging S{'!t, tt) and a non-zero value 
of Ps. For finite system sizes, both quantities are always 
finite and estimates for the thermodynamic limit are ob- 
tained by carefully studying the finite-size scaling of the 
observables. 

A jump in p with varying p indicates a discontinuous 
(first order) transition, and has been used to identify re- 
gions of phase separation in the canonical ensemble (fixed 
density p). We postpone a more rigorous analysis to ac- 
curately identify the nature of the transitions and the 
accurate domain boundaries to a later study, and focus, 
instead, on establishing the existence of SS phase over 
finite regions of the parameter space. 

A plot of p as a function of p shows clear indications 
of phase separation at p < 0.5 for all values of V — the 
discontinuity in p grows with increasing V. For p > 0.5, 
the curves are qualitatively different. For V < U/A, there 
is a small, but finite, region of positive slope (eg., 0.5 < 
p < 0.52 for V = 3), followed by phase separation for 
0.52 < p < 0.60, and a region of positive slope for p > 0.6. 
At V = U/A, there is no evidence of phase separation for 
p > 0.5. With V > U/4, the region of phase separation 
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FIG. 3: The average density as a function of the chemical po- 
tential for three different values of V {U — 20). For clarity of 
presentation, data for only one system size, L = 16, is shown. 
Error bars are smaller than the symbol sizes. Discontinuous 
transitions are marked by finite jumps in the particle density. 

shifts to large densities, p < 1. The location and extent 
of phase separated regions for small V{< U /A) agrees well 
with the results of Ref. pT| b apart from the extra region 
of positive slope for 0.5 < p < 0.52. As shown below, 
the ground state at these densities has supersolid order. 
The extent of the supersolid phase decreases rapidly with 
decreasing V, becoming vanishingly small in the limit 
V U. We note that for small V the excess density 
p — 1/2 < 1/L and larger L are required to map the SS 
phase boundary accurately. 

Ground state results for S{t:, tt) and ps as a function 
of p for three representative values of V are shown in 
Fig. ^ for three different system sizes L. The data are 
seen to be well converged with system size. At small p, 
the ground state is a superfluid (SF) — the stiffness con- 
verges to a finite value while ^(Tr, tt) scales to zero. As 
the density increases beyond a critical value Uci, there 
is a discontinuous transition to a (tt, tt) ordered charge- 
density-wave (CDW) ground state with p = 0.5. Any 
intermediate density is inaccessible in the grand canoni- 
cal ensemble. For V — 3t{< J7/4), at p > 0.5, there are 
indications for a small region of supersolid (SS) charac- 
terized by finite values of both Sijr, tt) and ps, but further 
finite size scaling tests will be needed to check whether 
this region remains in the thermodynamic limit. With 
increasing density, there is another discontinuous transi- 
tion to an SF state with a second region of PS. Finally, 
at p = 1, the ground state is a Mott insulator (MI) with 
both ^(Tr,?!) = and ps 0. For V = 5t(= C//4), 
the extent of the SS region increases substantially and 
its stability is well established. Additionally, the second 
phase separated region shrinks to zero and there is a di- 
rect SS-SF transition. For V = 6t{> C//4), the SF phase 
at high densities is replaced by another region of phase 
separation. The ground state at p = 1 changes from a 
MI to a (tt, tt) ordered CDW state with two particles oc- 
cupying every alternate lattice site, with a discontinuous 
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FIG. 4: The scaled static staggered structure factor (open 
symbols) and superfiuid stiffness (filled symbols) as a function 
of average particle density p for three representative values of 
V. At small p, the ground state is a superfluid (SF) for all V. 
With increasing density there is discontinuous transition to a 
(tt, 7r)-ordered charge density wave (CDW) at p = 0.5, with an 
intermediate region of phase separation (PS). At p > 0.5, the 
ground state is a supersolid (SS) with finite values of both 
^(Tr, tt) and ps- For V < (7/4, this is followed by another 
discontinuous transition to a SF (the precise mapping of the 
SS/PS boundary requires larger L), with a second region of 
PS. For V = [7/4, the SS region extends to higher densities 
with a continuous transition to a SF ground state at p « 
0.74. The SS state extends to even higher densities for V > 
(7/4, but is followed by a discontinuous transition to a (tt, tt)- 
ordered CDW at p = 1, with an accompanying region of PS. 
For y < (7/4, the ground state at p = 1 is a Mott insulator. 
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FIG. 5: The ground state phase diagram in the V — fi plane, 
notations as in Fig. |l|. Different sohd-ordered phases are 
shown schematically. The PS regions are manifested as dis- 
continuous transitions across the corresponding phase bound- 
aries (not shown). 

transition separating it from the SF state. 

The resuhs are combined to map the schematic ground 
state phase diagram of the Hamiltonian (|l|) in the V — 
p (Fig. |l|) and p. ~ 1/V (Fig. ||) planes. Fig. |l| shows 
the different phases in the (V, p) parameter space at a 
constant value of the on-site interaction, U — 20, t — 
1. For small V, the ground state is a SF for all p < 
1. At V > 2.5t, the different phases appear as shown 
in the figure. The extent of the supersolid phase and 
that of the phase separated region at p < 0.5 increases 
with increasing V, whereas the phase separated region 
at p > 0.5 gets vanishingly small for moderate values 
of V. It is not clear from the available data if the PS- 
SF phase boundary meets the SS boundary at a point, 
or approaches it asymptotically. At V > 5.0t, the SF 
region at high densities is replaced by a phase separated 
region, while the ground state at p — 1 changes from a 
MI to a (tt, 7r)-ordcred CDW with two bosons occupying 
every other lattice site. It should be emphasized that the 
phase diagram is qualitative and the phase boundaries 
are approximate. 

The features of the phase diagram as a function oil/V 
(Fig. ^ are markedly different from the "lobe" structure 
observed in a plot of /i as a function oit/U for the EBHM. 
The nature of the ground state a,t p — \ changes from a 
CDW to an MI as V is varied across U/A. This is ac- 
companied by a change in the curvatures of the phase 
boundaries. Furthermore, the MI region remains finite 
even in the limit oi V 0. No evidence of supersolid 
phase is found at p = 0.5, in agreement with the varia- 
tional studies and previous numerics |^ . 

In conclusion, we have used exact strong-coupling ex- 
pansion and QMC simulations to study the nature of 
the ground state phases of the extended boson-Hubbard 
model on a square lattice. The interplay of the on-site 
and nearest-neighbor interactions leads to a rich phase 



diagram including a supersolid phase with simultaneous 
diagonal and off-diagonal long-range order. We have 
provided strict arguments why a soft-core model with 
V > U / z and densities /? > 1/2 is sufficient to stabilize a 
supersolid phase in a model with nearest neighbor cou- 
plings only. This is in contrast to the "hard-core" bosons 
where the system phase separates for all values of the 
nn interaction strength and additional nnn interactions 
or hoppings are needed to stabilize a supersolid ground 
state. Also, in the studied range of parameters (includ- 
ing very large V < 12t) we have not found any nominally 
gapless phase with both ^(Tr, tt) and ps zero which could 
potentially be identified with a bose metaljo (see also 
the argument against such a phase in Ref . |ll|d) . 
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out on the computer cluster at the Institute of Geophysics 
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